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We analyse the number density and radial distribution of substructures and satel- 
lite galaxies using cosmological simulations that follow the gas dynamics of a baryonic 
component, including shock heating, radiative cooling and star formation within the 
hierarchical concordance ACDM model. We find that the dissipation of the baryons 
1^ ■ greatly enhances the survival of subhaloes, expecially in the galaxy core, resulting in 

| a radial distribution of satellite galaxies that closely follows the overall mass distribu- 

tion in the inner part of the halo. Hydrodynamical simulations are necessary to resolve 
\Q ' the adiabatic contraction and dense cores of galaxies, resulting in a total number of 

\ satellites a factor of two larger than found in pure dark matter simulation, in good 

•n ■ agreement with the observed spatial distribution of satellite galaxies within galaxies 

and clusters. Convergence tests show that the cored distribution found by previous 
authors in pure N-body simulations was due to physical overmerging of dark matter 
O ^ only structures. 

We proceed to use a ray-shooting technique in order to study the impact of these 
additional substructures on the number of violations of the cusp caustic magnification 
relation. We develop a new approach to try to disentangle the effect of substructures 
from the intrinsic discreteness of N-Body simulations. Even with the increased number 
of substructures in the centres of galaxies, we are not able to reproduce the observed 
high numbers of discrepancies observed in the flux ratios of multiply lensed quasars. 
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1 INTRODUCTION 

The study of substructures within simulated galaxy sized 
haloes have posed interesting problems to the current, 
widely accepted, ACDM scenario. As pointed out previ- 
ously, (Klypin et al 1999, Moore et al. 1999) the number 
of surviving subhaloes found in N-Body simulations greatly 
exceeds the number of observed satellites around the Milky 
Way and Andromeda. The properties of subhaloes on dif- 
ferent scales has been the subject of many recent studies 
which have pushed the resolution of dissipationless simula- 
tions (Moore et al. 1998, Ghigna et al. 1998, Colin et al. 
2000, Ghigna et al 2000, Springel et al. 2001, De Lucia et 
al. 2004, Kravtsov et al 2004, Diemand et al 2004a, Gao et 
al 2004a, Reed et al 2005, Zentner et al 2005a, Zentner et 
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al 2005b). The kinematical properties of subhaloes are now 
well understood - they make up a fraction of between 5 and 
10% of the mass of virialized haloes on scales relevant to 
observational cosmology. 

Most of these previous studies used dissipationless cos- 
mological simulations; although non-baryonic dark matter 
exceeds baryonic matter by a factor of Qdm/^b — 6 on aver- 
age, the gravitational field in the central region of galaxies 
is dominated by stars and gas. In the hierarchical galaxy 
formation model, stars are formed by the condensation of 
cooling baryons at the halo center. The cooling baryons in- 
crease the density in the central halo region mainly because 
of the extra mass associated with the inflow, but also be- 
cause of the adiabatic contraction of the total mass distri- 
bution (Blumenthal et al (1986), Loeb & Peebles (2003), 
Gao et al. (2004b), Gnedin et al. (2004)). This process acts 
both for the host halo and its subhaloes, therefore we may ex- 
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pect that the dark matter substructure haloes formed within 
hydrodynamical simulations will experience a different tidal 
force field and they themselves will be more robust to tidal 
effects. Since overmerging of dark matter substructures is 
sensitive to their central structure (Moore, Katz & Lake 
f996), pure N-Body simulations may bias results because 
of physical overmerging (Diemand, Moore & Stadel 2004a, 
DMS04 hereafter). 

The fact that substructure haloes are spatially more ex- 
tended than the averaged mass distribution was first pointed 
out by Ghigna et al. (f998). By increasing the resolution by 
over an order of magnitude, DMS04 showed that this result 
was independent of the numerical resolution. Furthermore 
the distribution of galaxies in clusters appears to follow the 
overall mass distribution, quite unlike the subhaloes selected 
by final bound mass in pure dark matter simulations. By 
tracing haloes backwards and forwards in time through the 
simulations they argued that the missing central subhaloes 
may be destroyed by tidal forces at redshifts higher than 
z—5. The survival of substructure and galaxies within dense 
environments has implications for indirect detection tech- 
niques, for example using image distortions of gravitation- 
ally lensed distant quasars by foreground galaxies and their 
dark matter haloes. 

ft has been argued that a possible signature of the pres- 
ence of dark matter substructures can be found in strong 
gravitational lensing of QSOs (Mao & Schneider f998; Met- 
calf & Madau 200f; Chiba 2002; Metcalf & Zhao 2002; 
Dalai & Kochanek 2004; Chen, Kravtsov & Keeton 2003; 
Kochanek & Dalai 2004; Mao et al 2004, Amara 2004). if 
a distant image source is close to (or inside) a cusp in a 
caustic curve, three of the images will be clustered together 
and the sum of their magnifications will be zero (Zakharov 
(f995), taking the negative parity image to have negative 
magnification). This relation holds for a wide class of smooth 
analytic lens models (i.e. Schneider & Weiss f992, Keeton 
et al. 2003); on the other hand all known observed lensed 
QSOs violate this relation which has been explained as due 
to the presence of cold dark matter substructure within the 
lensing galaxy's halo. However, some of these discrepant sys- 
tems may be due to microlensed stars rather than to cold 
dark matter substructure (Keeton et al 2003). Bradac et al. 
(2004) analysed a low resolution simulation of a galaxy and 
claim that the level of substructure present in simulations 
produces violations of the cusp relation comparable to those 
observed. Amara et al. (2004) implanted an idealised model 
of galaxy into the center of a high resolution galactic halo 
extracted from dissipationless N-Body simulations to test 
the effects of substructure on lensed images. Their findings 
contrast those of Bradac et al 2004, since they found that the 
substructures produced in a ACDM halo are not abundant 
enough to account for the cusp caustic violation observed, 
these results are also confirmed by Mao et al 2004, based 
on analysis of substructure abundance in pure dark matter 
simulations. 

The first part of this paper is devoted to the analysis 
of the subhalo population around a galaxy sized halo that 
forms in a large cosmological simulation simulated with dark 
matter only and then with the inclusion of a baryonic com- 
ponent. In section 2 we present the numerical simulations, 
our halo-finding scheme and resolution tests. The proper- 
ties of the main halo and its satellites are presented in sec- 



tion 3. Comparison with observations and a discussion of 
the results are in section 4. In the second part of this pa- 
per we re-examine the effects of substructures on multiply 
lensed quasar images using our new hydrodynamic galaxy 
simulation in combination with a ray shooting technique. In 
section 5.1 we present the lensing code, while section 5.2 is 
devoted to multiple images analysis and in sections 6 and 7 
we present our results and the conclusions of our work. 



2 NUMERICAL SIMULATIONS 

The simulations were performed with GASOLINE, a multi- 
stepping, parallel TreeSPH A-body code (Stadel 2001, Wad- 
sley et al. 2004). We include radiative and Compton cooling 
for a primordial mixture of hydrogen and helium. The star 
formation algorithm is based on a Jeans instability criteria 
(Katz 1992), where gas particles in dense, unstable regions 
and in convergent flows spawn star particles at a rate pro- 
portional to the local dynamical time (see also Governato et 
al 2004). The star formation efficiency was set to 0.1, but in 
the adopted scheme its precise value has only a minor effect 
on the star formation rate (Katz 1992). The code also in- 
cludes supernova feedback as described by (Katz 1992), and 
a UV background following Haardt & Madau 1996. 

We have selected a candidate Galactic mass halo 
(Md m ~ 10 12 Mq) from an existing low resolution dark 
matter only simulation in a concordance (A=0.7,S7o=0.3, 
(78=0.9) cosmology and resimulated it at higher resolu- 
tion using the volume renormalization technique (Katz & 
White 1993), and including a gaseous component within 
the entire high resolution region. The mass per particle 
of the dark matter and gaseous particles are respectively 
m d = 1.66 x 10 6 M Q and m g = 3.28 x 10 5 M Q . The dark 
matter has a spline gravitational softening length of 200 pc 
and we have about 1 x 10 6 particles for each component 
(dark and gas) in the high resolution region. 

The SPH simulations with cooling are computationally 
expensive and we are forced to stop the full calculation once 
the galaxy has formed, at a redshift z — 1.5. (The parallel 
calculation is dominated by the few remaining gas particles 
which need extremely small timesteps in order to satisfy the 
Courant criterion). In order to study the dynamical evolu- 
tion of the stellar and dark matter satellites we have evolved 
the simulation to the present epoch without following the re- 
maining gaseous particles - we turn them into collisionless 
particles and treat only their gravitational interactions. We 
do not believe that this influences any of our conclusions 
since (i) most of the gas has already turned into stars be- 
tween a redshift z=5.5 and z=2.5, and (ii) continuing to 
include cooling of the remaining gas would only allow us to 
resolve higher density gas clouds. Inside the virial radius at 
a redshift z=0 we have 8.6 x 10 5 dark matter particles and 
2.4 x 10 5 star particles. In Figure 1 we show the star forma- 
tion rate in our simulation which is in good agreement with 
previous numerical studies (i.e. Governato et al. 2004). 

The same object has also been simulated with dark mat- 
ter only using the same spatial and mass resolution as the 
hydro run, in addition one simulation with four times bet- 
ter mass resolution (this object is Gall in DMS04). In the 
following we will refer to the three simulations as hydro, dm 
and drriHR respectively. 
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Figure 1. The star formation rate in the hydro simulation. The 
time since the big-bang is indicated and the red vertical line cor- 
responds to 2 = 2, where cooling were switched off. 

3 HALOES PROPERTIES 

3.1 Host halo 

In Figure 2 we show the rotation curve (defined as V c (r) = 
\J GM{< r)/r) for all components, stars and dark matter 
separately for the galaxy at z=0. For comparison we also 
show the dark matter rotation curve obtained in the pure 
dark matter simulation of the same object. The baryons 
dominate within the inner 10 kpc and the effect of the adi- 
abatic contraction on the dark matter can be clearly seen 
- the baryons have increased the mass within 5 kpc by a 
factor of four. The steeper central cusp can be clearly seen 
in the density profile plot (fig. 3), the dm halo has a profile 
that can be well fitted by an NFW profile with a concentra- 
tion c v ir = 9.6 (defined with respect to the virial radius); on 
the other hand in the presence of gas and stars the density 
profile in the inner region has an almost constant cusp slope 
q « —2.0, as predicted by the adiabatic contraction model 
(Blumenthal et al. 1986). In figure 1 we show the star for- 
mation rate (SFR), the bulge of the galaxy forms through a 
series of rapid major merger events that end around z — 2.8, 
turning mostly low angular momentum gas into stars. There 
are no more merging events in the formation of our galaxy 
after z = 2 so even if we cannot follow the SFR directly 
beyond z = 1.5 we do not expect it to be different from a 
slowly decreasing function of time. 

3.2 Sub haloes 

Within the virial radius of the high resolution CDM simula- 
tions we can resolve several hundreds of substructure haloes 
(bound over-dense clusters of particles). We identify sub- 
haloes with SKID (Stadel 2001), which calculates local den- 
sities using an SPH kernel, then moves particles along the 
density gradient until they oscillate around a point (i.e. move 
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Figure 2. Rotation curve (defined as V C (R) = y/GM(< R)/R) 
for all components (upper solid line), star (dotted line) and dark 
matter (dashed line) for the galaxy at z=0. For comparison the 
dark matter rotation curve (lower solid line) obtained in the pure 
dark matter simulation of the same object is also shown. 



less than some length I). Then they are linked together using 
FoF with this I as a linking length. 

SKID with / = 4eo (where to is the gravitational soft- 
ening of the simulation) adequately identifies the smallest 
subhaloes and the centers of the largest subhaloes. For the 
latter the calculated bound mass is underestimated. Using 
/ = 10eo can cure this, but then some of the small subhaloes 
are missed. Therefore we have used a combination of the 
subhalo catalogues obtained with these two linking lengths 
in order to create the complete catalogue of subhaloes and to 
calculate their correct structural parameters. We included in 
our study only halos with M > 2 x 1O 8 M (i.e. ridm > 100). 

In Figure 4, where we plot the dark matter subhaloes 
mass function in the hydro and dm simulations, calculated 
within a sphere of radius 1 Mpc centred on the galaxy. As 
expected, we do not observed a big difference between the 
two simulations in such a large volume, this is because the 
presence of baryons is thought to be much more important 
in the inner region of our galaxy where the tidal field is much 
more efficient in destroying subhaloes. The effect of baryons 
can be easily seen in Figure 5, where we plot the ratio of 
the number of subhaloes in the two simulations as function 
of the distance from the center of the galaxy. At 70 kpc (~ 
1 /3i?200 ) we see an increase by a factor of two in the numbers 
of surviving satellites. Within the inner 40 kpc the number 
of satellites with mass greater than 2 x 1O 8 M0 is enhanced 
by a factor 4-5 over that found in the pure dark matter 
simulation. We found that this difference is more striking 
for the most massive subhaloes, the number of satellites with 
mass greater than 10 7 M© does not differ so much from the 
N-body result. This may be revealing the limitation of our 
resolution - the smallest haloes of mass ~ 10 7 M Q are just 
resolved with about 10 DM particles and are more easily 
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Figure 3. Dark matter density profile in the N-Body (dashed 
curve) and Hydro (solid line) simulations, the effect of the adi- 
abatic contraction can be clearly seen in the inner part of the 
profile. 

destroyed through numerical effects. Another explanation 
could be the inhibition of star formation in smaller haloes 
because of the ultra-violet (UV) background and supernovae 
(SN) feedback that raise the temperature of the gas and stop 
its collapse in the DM halo if it is not sufficiently bound 
(i.e. its potential well is not deep enough). Without stars, 
the DM can be more easily destroyed by the tidal forces of 
the main halo (see also section 4 for more details on the 
feedback effects). 

As noticed by previous authors (DMS04 for a recent 
analysis) the spatial distribution of subhaloes in cold dark 
matter simulations of galaxies is anti-biased with respect to 
the mass; Nagai & Kravtsov (2005) have recently shown that 
part of this bias is due to the varying amount of mass loss at 
different radii, and that it is considerably smaller if instead 
of using the mass of satellites at z — 0, one uses the mass 
measured at the accretion time. For our purposes we have 
decide to use the mass at the present time, because we are 
interested in the different tidal forces and mass loss between 
the hydro and dm run. 

Figure 6 shows the relative number density of surviving 
satellites M > 2 x 10 8 at z=0, in several pure dark matter 
simulations and our hydro simulation. In all dm cases, the 
satellites are more extended than the overall mass distribu- 
tion (dominated by the smooth dark matter background). 
This effect becomes larger inside of half the virial radius. 
Within this region, the dark matter only simulations reveal 
a much natter number density distribution. The satellite dis- 
tribution in the hydro simulation is much steeper and it is 
rather similar to the smooth mass distribution even if also 
in this case the overall radial distribution of the subhalos 
is more extended than the underlying dark matter distribu- 
tion as shown clearly by the excess in the number density of 
subhalos at 0.5 < r/r V i r < 1.0. 
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Figure 4. Dark matter subhaloes mass function in the hydro 
(black) and dm (blue) simulations, in a sphere of 1 Mpc around 
the center of the galaxy. 

The inner core distribution found in the dm run almost 
disappears and the satellite profiles are well fitted by an 
NFW like function even if they are still less concentrated 
than the overall mass distribution (in agreement with find- 
ings from hydro simulations on cluster scales by Nagai & 
Kravtsov (2005)); the concentration parameter (r V i r /r c ) for 
the satellite distribution is ~ 6.5, where the one for the 
smooth DM distribution in the range 0.07 < r/r V i r < 1.0 is 
cdm = 9.6. 

3.3 Resolution tests 

In the hydro simulation we have a higher mass resolution 
due to the presence of baryons; the mass of the star and gas 
particles is roughly 20% of the mass of the DM particles, this 
means that in principle the lower number of subhaloes close 
to the center that we found in the N-Body simulation can be 
partially due to an overmerging phenomenom that depends 
on the resolution. In the lower panel of Figure 7 we compare 
the number of subhaloes as a function of radius between 
the hydro simulation and the highest resolution dark matter 
only run (drriBR) which has a mass resolution a factor 4 
better than the low resolution case. Figure 7 demonstrates 
that the increase in the number of substructure haloes in 
the hydro run is apparent even when compared to the higher 
resolution N-Body simulation. 

The second test that we have performed is related to 
the force resolution (gravitational softening) adopted for the 
star particles. Stars are concentrated in the inner few kpc 
of the galaxy and a small softening must be chosen to cor- 
rectly follow the dynamical evolution of such high density 
regions. We have run the hydro simulation for two different 
values of the star gravitational softening (e): ei = 0.2 kpc 
and £2 = 1.5 kpc. Results are presented in the upper panel 
of Figure 7 where we show the difference in the number of 
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Figure 5. Ratio between the number of DM subhaloes with M > 
2 X 10 8 Mq in the hydro and dm runs as function of the distance 
from the center of the main halo. 
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Figure 6. Substructure radial density profiles for different simu- 
lations. 

subhaloes as a function of the distance from the center. The 
effect of softening is important - too large a value leads to 
tidal disruption of satellite galaxies since they have artifi- 
cially shallow central potentials and are hence less bound. 
The larger stellar softening (e2 = 1.5 kpc) erases the stabi- 
lizing effect of the stars inside the satellites and the radial 
distribution of surviving subhaloes is identical to the pure 
dark matter case (see the dashed line in the upper panel of 
Figure 7). 



Figure 7. Effect of resolution on the number of subhaloes. Up- 
per panel: ratio between the number of DM subhaloes (M > 
2 X 1O 8 M0) in the hydro run for two different values of the star 
particles softening: ei = 0.2 kpc and t2 = 1-5 kpc as function of 
distance from the center. The dashed line is the same quantity 
but for the hydro run with e = £2 and the dm run. Lower Panel: 
the ratio between the number of DM subhaloes (M > 1O 7 M0) 
in the hydro run and in the high resolution dm^R (dashed line) 
and low resolution dm (solid line) runs. 

4 COMPARISON WITH OBSERVATIONS 

For this comparison we used the data for the Local Group 
as compiled by Mateo et al. (1998). In this sample only rela- 
tively massive satellites with an estimated rotational veloc- 
ity or a three-dimensional velocity dispersion of stars greater 
than 10 km sec" 1 are considered. In order to simplify the 
comparison even further, as Klypin et al. (1999) we have 
considered the number of satellites within a radius of 280 
kpc from the Milky Way and Andromeda, which is close to 
the expected virial radius of these galaxies. 

Even though the remaining gas particles are turned into 
collisionless tracers at z = 1.5, this does not affect the kine- 
matical properties of the satellites at z=0. All of the satel- 
lites that end up within the virial radius of the galaxy halo 
have formed prior to this epoch. Subsequent loss of the re- 
maining gas particles to stars or by stripping processes would 
not lead to significant changes to the satellite structural pa- 
rameters. 

In Figure 8 we show the cumulative velocity distribution 
function (VDF) of satellites: the number of satellites per 
unit volume and per central object with internal circular 
velocity larger than a given value V c . We show the VDF for 
different runs: the solid line represents the results for the dm 
simulation, as already known the CDM models over-predict 
the number of satellites for V c < 30 km sec -1 . The dashed 
line represent the results from Klypin et al 1999, that were 
originally compared with the data from Mateo et al 1998. 
The simulation results agree quite well above 20 km sec -1 , 
our higher resolution pure dark matter simulations allow 
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Figure 8. Cumulative circular velocity distribution function of 
satellites. Black solid line represents the results for the DM only 
simulation, blue solid line is the VDF for stars (visible) haloes, 
the dashed blue line is for DM haloes. Black triangles with error 
bars show average results for Milky Way and Andromeda Satel- 
lites (Mateo 1998 & Klypin et al 1999). The red dashed line are 
results for a ACDM dissipationless galaxy obtained by Klypin et 
al (1999). 



us to follow the distribution of subhaloes to lower circular 
velocities (V c > 10 km sec -1 ). 

The open circles connected by a solid line shows the 
VDF of dark matter haloes for the hydro run. The numbers 
of satellites in the range 15 < V c < 30 km sec -1 is in- 
creased, and now also for V c > 35 km sec -1 there is an over 
abundance of satellites. We have also computed the VDF 
for satellites in the hydro run considering only their stars, 
shown by squares connected with the long dashed line. This 
is a fairer comparison since we are comparing stellar kine- 
matics in each case. Stars are more concentrated than DM, 
so they trace the dynamics within a smaller central region, 
this explains why the VDF for stars is below the DM one 
over the whole velocity range (Hayashi et al 2004) but this 
effect is much to small to reconcile the simulations with the 
observations of Local Group satellites. 

Another clear difference between the two VDFs is that 
we have no stellar satellites with V c < 15 km sec -1 , which 
is due to the supernovae feedback and UV background in 
the simulation. To confirm this we have run a full hydrody- 
namical simulation without these external feedback sources, 
shown in fig. 9, where the VDF for the no feedback case 
continues to rise within V c < 15. With the weak feedback 
used in this paper hydrodynamic cosmological simulations 
produce VDF's which lie above the pure dark matter results 
and the discrepancy to Local Group observations becomes 
even larger. Runs with strong feedback from reionization are 
able to produce realistic VDF's, we will present such runs 
in a forthcoming paper (Maccio et al. 2005b) . 
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Figure 9. A similar plot to figure 8: the solid (blue) line is the 
stellar satellite VDF in absence of SN feed back, the dashed line is 
with feedback turned on, black dotted line represents the results 
for the DM only simulation. 

5 LENSING ANALYSIS 

Quasars that are being gravitationally lensed into multiple 
images have recently been used to place limits on the sur- 
face density of cold dark matter subhaloes (Mao & Schnei- 
der 1998, Metcalf & Madau 2001, Chiba 2002, Metcalf & 
Zhao 2002, Dalai & Kochanek 2002, Chen Kravtsov & Kee- 
ton 2003, Bradac et al. 2004, Mao et al 2004, Amara et al. 
2004) . Small mass clumps that happen to lie near the images 
affect the observed magnification ratios. The question arises 
as to whether these observations are compatible with distor- 
tions expected to occur from dark matter substructures and 
satellite galaxies within the ACDM model. In the following 
we will first present our lensing code, then we will briefly re- 
call the main features of the so called "Cusp Relation" and 
finally we will report our results. 

5.1 Lensing Simulations 

Our ray shooting code is described in detail in Maccio 2005, 
here we will only summarize its mean features. The galaxy 
is centered in a cube of length 0.6 Mpc and we study three 
lense images, obtained by projecting the particle positions 
along the three coordinate axes. We then divide the pro- 
jected density field E by the critical surface mass density 
for lensing 

_ c 2 D S 

thus obtaining the convergence k. Here c is the speed of 
light, G is the gravitational constant, while Dl, Ds, Dls are 
the angular-diameter distances between lens and observer, 
source and observer, lens and source, respectively. In the 
following, we adopt zl = 0.3 for the lens redshift and z s = 
3.0 for the source redshift. 
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The deflection angle due to this 2D particle distribution, 
on a given point x on the lens plane reads: 



N 

E 

2=1 



4G 



(2) 



Here yj is the position of the j-th particles and N is the 
total number of particles. 

Since a direct summation requires a long time, we speed 
up the code by using a P 3 M-like algorithm: the lens plane 
was divided into 256x256 cells and direct summation was 
applied to particles belonging to the same cell of x and for 
its 8 neighbor cells. Particles in other cells were then seen as 
a single particle in the cell baricenter, given the total mass 
of the particles inside the cell. The deflection angle diverges 
when the distance between a light ray and a particle is zero. 
To avoid this unwanted feature we introduce a softening 
parameter, e in equation (2); the value e is tuned to the 
resolution of the current simulation. 

We start to compute a(x) on a regular grid of 
4096x4096 test rays that cover the central quarter of the 
lens plane (multiple images form very close to the lens cen- 
ter). This resolution does not provide enough pixels in the 
inner regions affected by strong lensing to model lensing 
properties in the correct way. The resolution is increased by 
extracting the central region where strong lensing is occur- 
ring and using bilinear interpolation to calculate the rele- 
vant quantities to higher resolution. Our final resolution is 
equivalent to a bundle of 16384x16384 light rays. 

The relation between image and sources positions is 
given by the lens equation: 



y 



a(x) 



(3) 



and the local properties of the lens mapping are then de- 
scribed by the Jacobian matrix of the lens equation, 

da h 
dxk 

and the magnifications factor [i is given by the Jacobian 
determinant of A, 

1 



A (#\ Q V h X 

Ahk{x) = - — = dhk 

OXk 



(4) 



[A 11 (x)A 2 2(x) - A 12 (x)A 2 i(x)]~ 



(5) 



det A 

Finally, the Jacobian also determines the location of 
the critical curves x c on the lens plane, which are defined 
by det A(x c ) = 0. Because of the finite grid resolution, we 
can only approximately locate them by looking for pairs of 
adjacent cells with opposite signs of det A Then the lens 
equations 



y c — x c — a(x c ), 



(6) 



yields the corresponding caustics y c , on the source plane. 

To find the images of an extended source, all image- 
plane positions x are checked to see if the corresponding 
entry in the map table y lies within the source: i.e. for a 
circular source with radius r c and centered in (j/f; y 2 ) it is 
checked if: 



■Vi 



(2/2 



e\2 2 

V2) <r c , 



(7) 



where (yi,y 2 ) are the components of the vector y. The 
sources are modeled as circles with a radius of 60 pc ac- 
cording to Amara et al (2004). Those points fulfilling the 
previous equation are part of one of the source images and 
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Figure 10. Basic lens configuration. The caustic surface is shown 
as a black line and critical curves are shown as cyan lines. The four 
images that are usually observed are shown as red dots, and the 
three that are used to inspect the cusp relations are the ones inside 
the small black circles. The blue square is the source position. The 
softening adopted for this lens map is 0.5 Kpc. 



are called image points. We then use a standard friends- 
of-friends algorithm to group together image points within 
connected regions, since they belong to the same image. 

A typical lens configuration is shown in Figure 10 where 
critical lines, caustic lines, source and images position are 
indicated. The images within the black circles correspond to 
the those selected for the cusp relation investigation. 



5.2 Cusp Caustic Relation 

There are basically three configurations of four-image sys- 
tems: fold, cusp, and cross. In this paper we will mainly 
concentrate on the cusp configuration, that corresponds to 
a source located close to the cusp of the inner caustic curve. 
The behavior of gravitational lens mapping near a cusp was 
first studied by Blandford & Narayan (1986), Schneider & 
Weiss (1992) and Mao (1992) and Zakharov (1995), who in- 
vestigated the magnification properties of the cusp images 
and concluded that the sum of the signed magnification fac- 
tors of the three merging images approaches zero as the 
source moves towards the cusp. In other words (e.g. Za- 
kharov 1995) : 



Reus 



HA + [IB + HC 

\^a \ + + \(ic\ 



0, for Htot 



(8) 



where [it ot is the unsigned sum of magnifications of all four 
images, and A,B & C are the triplet of images forming the 
smallest opening angle. By opening angle, we mean the angle 
measured from the galaxy center and being spanned by two 
images of equal parity. The third images lies inside such 
an angle. This relation is an asymptotic relation and holds 



© 2005 RAS, MNRAS 000, 000-000 



8 A.V. Maccid, B. Moore, J. Stadel & J. Diemand 



R ratio 




arc sec 

Figure 11. Value of the quantity R C usp in the four images region 
of the source plane, here a softening of 0.3 Kpc is used. 

when the source approaches the cusp from inside the inner 
astroid caustic. 

Since we know the lens position and the source posi- 
tion the procedure of finding the cusp images is straightfor- 
ward, we have identified the triplet of images belonging to 
the smallest opening angle, we have seen that the cusp im- 
ages are better identified using as opening angle measured 
from the fourth image and being spanned by two images of 
equal parity. We have used this method only to find the cusp 
images, instead for testing the cusp relation (see eq. 9) we 
have used the opening angle as defined from the center of 
the galaxy (A9). 

Approximately 25000 lens systems are generated with 
the source position inside the astroid caustic. Figure 11 in- 
dicates in color the value of R cusp for all the different source 
positions; the apparent discontinuities originate from differ- 
ent image identifications. In the very center of the caustic 
the cusp relation is not well defined (what you have is mostly 
a "cross relation", four images situated at the vertices of a 
cross centered on the cusp center), as the source moves in 
the direction of the minor or major axes we choose different 
subsets of three cusp images and therefore the discontinuity 
arises. 



6 EFFECTS OF SUBSTRUCTURES 

Because of the finite size (discreteness) of the particles in 
the simulations there is a significant amount of shot noise in 
the surface density estimate, which can affects the lensing 
properties. The usual approach (Bradac et al. 2004, Amara 
et al. 2004) is to use a Gaussian kernel to smooth the surface 
density, in Amara et al. (2004) a detailed study of the im- 
pact and calibration of the smoothing was presented, their 
main conclusion was that a smoothing of 0.5 kpc is suitable 
for studying the cusp relation properties without losing any 
spatial information. 



In order to try to disentangle the effects of bound sub- 
structures to the spurious effect produced by single finite 
particles, we have adopted a novel approach. We have have 
tried to remove substructures from the lensing halo to see 
how this can change the cusp relation. Our removal pro- 
cedure works as follows: first we identified all bound sub- 
structure using SKID (see section 3.2), then each particle 
belonging to any of the subhaloes is rotated around the cen- 
ter of the galaxy using three random Euler angles. In order 
to avoid thin circular shells of particles we added a random 
±10% error to each distance. We want to emphasize that 
we do not physically remove any substructures, because this 
will change the overall properties of the lens: mass, density 
profile, etc. 

This procedure allows us to smooth out only the sub- 
structures leaving unaltered all the main features of the pri- 
mary lens. In Figure 14 (upper panel) we show the integral 
radial mass profile before and after the removal of the sub- 
structures, not only the total mass is preserved but also its 
radial profile, this means that all the lensing properties of 
the galaxy remain unchanged (critical curves, caustics, posi- 
tion of images) allowing us to make a one to one comparison 
between the two lensing maps. In Figure 12 and Figure 13 
the density map of the galaxy is shown, with and without 
substructures respectively. 

The cusp relation defined by equation 8 holds when the 
source is close to the cusp. As soon as the source moves away 
from the cusp deviations from R C us P = are observed, even 
for the smooth lens model. On the other hand the closer the 
source is to the cusp, the smaller is the angle spanned from 
the three images, so in order to take into account also the 
position of the source in evaluating the cusp relation it is 
better to define Rc Uap as (see also Amara et al. 2004): 

Rcusp -^jRcusp- (9) 

where A8 is the opening angle spanned by the two images 
with positive parity defined from the center of the galaxy. 

With this new definition of Rc Usp a set of three images 
is said to violate the cusp relation if Rc Usp > 1. This makes 
the comparison between simulations and observations much 
more straightforward. 

The differences in the reduced cusp relation violation 
in the two cases are shown in Figure 16, where we plot the 
number of sources that violate the reduced cusp relation as 
a function of the Gaussian smoothing scale (e s ). For both 
the lens models, the number of sources that violate equation 
9 decreases with the smoothing length because the effect of 
smoothing is to both reduce the impact of substructures and 
the noise introduced by single particles. 

The difference between the two is not so large, with a 
maximum for ec ~ 0.5 kpc, where the number of violations 
grows from 19% to 23%, this is because this value of eo is 
large enough to cancel the Shot noise, but not large enough 
to smear out the subhaloes in the simulation, in good agree- 
ment with the results of Amara et al. For smaller value of 
€g the signal is almost completely dominated by Shot noise, 
and for larger values we smooth too much, losing spatial 
information on the surface density of the lens. 

Figure 16 clearly shows that the impact of substruc- 
ture in a mass range 10 7 — 10 9 is very weak in disturbing 
the cusp relation. The resolution achievable with current 
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numerical simulations is still too poor to extend the analy- 
sis to a lower mass range of subhaloes, analytic arguments 
or semi-analytic prescriptions must be used (Maccio & Mi- 
randa 2005). 

Nevertheless a tentative comparison with observation 
can be made; there are 5 observed cusp caustic lenses sys- 
tems: B0712+472 (Jackson et al 1998), B2045+265 (Koop- 
mans et al 2003), B1422+231 (Patnik & Narasimha 2001), 
RXJ1131-1231 (Sluse et al. 2003) and RXJ0911+0551 (Kee- 
ton et al 2003); the first three are observed in the radio band, 
the last two in optical and IR. Three of them violate the re- 
duced cusp relation (i.e. R C us P > A9/2tt). This means a 
60% violation, that is significantly larger than the 15-25% 
we found for our simulations. The sources size used in this 
work (60pc) allows us to make a direct comparison mainly 
with QSO observed in radio than in optical or IR, even in 
this case we have a violation of the reduced cusp relation in 
2 objects over 3, that means 66% violation. Figure 15 shows 
the distribution of the values of the reduce cusp relation 
Rcusp both for data and simulates systems (with ec = 0.5 
kpc). Simulation results are unable to reproduce the high 
value tail that arises in the observational data; again it is 
possible to see that the effect of subhaloes is very weak in 
disturbing the cusp relation, and they only marginally en- 
hance the number of systems with R® uap > 1. 

As already discussed above, the difference between data 
and simulation results can arise from the current resolution 
limitations in the Numerical simulations, we expect to have 
many more small DM haloes close to the center of the galaxy 
with masses in the range 10 -3 — 10 6 M©, these haloes will be 
very concentrated and so they can more easily survive the 
tidal force of the central halo and if one of these haloes is 
close enough to one of the images it can perturb its magni- 
fication and so violate the cusp relation, (the smaller is the 
mass of the subhalo the closer it must be to the projected 
image position). 

Another possible explanation for the observed cusp re- 
lation violation, can be ascribed to the effect of the all sub- 
haloes that are in the galactic space along the line of sight 
of the lens (Metcalf 2004); the importance of these inter- 
galactic haloes with mass < 10 s Mq depends on the radial 
profile of the dark matter haloes and the primordial power 
spectrum at small scales. 

What is clear from our analysis is that the cusp relation 
violation can not be due to substructures in the primary lens 
with a masses above 10 7 Mq. This is true even though the 
inclusion of baryons has increased the projected numbers of 
subhaloes by a large factor. 



7 DISCUSSION AND CONCLUSIONS 

Using a high resolution hydrodynamic galaxy formation sim- 
ulation we have studied the number density and the spatial 
distribution of subhaloes within a Milky Way sized CDM 
halo. Baryons can concentrate at the halo centre via two 
mechanisms: The shock heated gas can radiate away energy 
causing the particles to fall inwards. Cold gas clouds can 
accrete along filaments, colliding and coalescing with ex- 
isting gas at the halo centre. The slowly infalling baryons 
has the effect of increasing the density of the dark mat- 
ter via adiabatic contraction. The overall effect of including 




Figure 12. Density map of the mass distribution within the full 
hydrodynamical simulation. The size of the box is 200 kpc. 




Figure 13. Density map of the smoothed mass distribution after 
randomizing the positions of the particles within the substructure 
haloes. 

baryons is to steepen the mass profile to nearly isothermal 
with p{r) oc r~ 2 from a shallower r _1 cusp. 

As already predicted by several (semi) analytic stud- 
ies, our simulations show that this central concentration 
of baryons enables subhaloes to better withstand the tidal 
forces generated by the main halo. This leads to an increase 
in the total number of subhaloes, especially within the inner 
third of the virial radius, such that they follow the overall 
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Figure 14. Mass versus radius before and after the removal of 
substructures, the overlap between the two curves is very good, 
in the lower panel the residuals are shown 



Figure 16. Fraction of total number of sources that violate the 
cusp relation as a function of the Gaussian smoothing length e 9 . 
Solid line is for the whole galaxy, dashed line for the galaxy with- 
out substructures (see text for its definition). 
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Figure 15. Distribution of R^ UBp values. The solid and dotted 
lines show the simulation results before and after the substructure 
removal. The long dashed line represents the observational data. 



mass distribution. This is in excellent agreement with the 
distribution of galaxies within galaxy clusters. However on 
the scale of Galactic mass haloes, this increases the discrep- 
ancy between the numbers of surviving satellites in CDM 
models and the flat luminosity function observed within the 
Local Group. A possible solution to this problem may come 
from the early reionization of the universe from early struc- 



ture formation (Madau & Rees 2001). Reionization raises 
the entropy of the gas that is required to fuel galaxy for- 
mation, preventing it from accreting into small dark matter 
haloes and lengthening the cooling time of that gas which is 
accreted (Bullock et al 2001, Somerville et al 2002, Ricotti, 
Gnedin & Shull 2002). Moreover as suggested by Benson & 
Madau (2003) winds from pre-galactic starbursts and mini- 
quasars may pollute the intergalactic medium (IGM) with 
metals and raise its temperature to a much higher level than 
expected from photoionization and so inhibit the formation 
of early galaxies. 

A detailed study of the impact of reionization on the 
formation of dwarf galaxies in hydrodynamic simulation will 
be presented in a forthcoming paper (Maccio et al 2005b), 
where we show that with an appropriate choice of the reion- 
ization parameters and modeling it is possible to reconcile 
observational data with the steep mass function of haloes 
and subhaloes within the cold dark matter model. 

In the second part of this work we have explored the 
consequences of the increased numbers of satellites for multi- 
ply lensed quasar images by foreground galactic mass haloes. 
In particular, the signatures on the violation of the so called 
cusp relation. Previous work have reached different conclu- 
sions on this issue: Bradac et al (2004), found an agreement 
between simulations and observations but their results were 
limited by the low resolution of their simulations. Amara 
et al (2004) (see also Mao et al 2004) have shown that is 
hard to reconcile the observed high number of cusp rela- 
tion violation with results from simulations moreover they 
have also show that Shot noise due to the discreteness of the 
simulation (where every particle is in principle a substruc- 
ture) plays an important role in changing the properties of 
the lens map. In order to disentangle the effect of substruc- 
tures from the effect of having an intrinsic discrete distribu- 
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tion of matter we compared results between haloes with and 
without any substructure present. We have developed a new 
technique to remove the satellites from the simulation with 
out changing the overall matter distribution of the primary 
lens. This analysis demonstrates that the impact on lensing 
of subhaloes in the mass range 10 7 - 10 10 Mq is very small. 
Also having a number of subhaloes which is about 8 times 
higher of the observed one in this mass range, the number 
of multiple lensed QSO that show a violation of the cusp 
relation (defined as eq. 9) is less than 24%, in contrast with 
an observed one of about 60%. Our results extend down to 
subhaloes masses of 10 7 M Q due to the resolution limit of our 
simulations. Even considering the impact of smaller masses 
for subhaloes using an analytic approach does not help in 
solving the cusp relation problem (Maccio & Miranda 2005) . 
A possible explanation could be that more variables must 
be taken into account in the lensing analysis, such as all the 
(sub)haloes that lie along the line of sight between us and 
the lens (see Chen et al. 2003, Metcalf 2004). 



ACKNOWLEDGMENTS 

The authors acknowledge J. Wadsley for development of the 
GASOLINE code and thank him for its use in this work and 
thank A. Kravtsov, D. Nagai, C. Mastropietro, L. Mayer, 
M. Miranda and P. Saha for useful discussions during the 
preparation of this work. We also thank the referee Hong- 
sheng Zhao for his comments. AM also thanks S. Kazantzidis 
for the use of the extract software. JD is supported by 
the Swiss National Science Foundation. All the numeri- 
cal simulations were performed on the zBox supercomputer 
(http://www-theorie.physik.unizh.ch/~stadel/) at the Uni- 
versity of Zurich. 

REFERENCES 

Amara, A, Metcalf, R.B., Cox, T.J. & Ostriker, J.P., 2004, astro- 
ph/0411587 

Blandford R., Narayan R., 1986, ApJ, 310, 568 
Bertschinger E., 2001, ApJSS, 137, 1 
Benson A. J., Madau P., 2003, MNRAS, 344, 835 
Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, 
ApJ, 301, 27 

Bradac M., Schneider P., Lombardi M., Stcinmctz M., Koopmans 
L. V. E., Navarro J. F., 2004, A&A, 423, 797 

Bullock J. S., Kravtsov A. V., Weinberg D. H., 2000, ApJ, 539, 
517 

Chen J., Kravtsov A. V., Kccton C. R., 2003, ApJ, 592, 24 
Chiba, M. 2002, ApJ, 565, 17 

Colin P., Klypin A. A., Kravtsov A. V., 2000, ApJ, 539, 561 

Dalai, N. & Kochanek, C. S. 2002, ApJ, 572, 25 

De Lucia C, Kauffmann C, Springel V., White S. D. M., Lanzoni 

B., Stoehr F., Tormen G., YoshidaN., 2004, MNRAS, 348, 333 
Diemand J., Moore B., Stadel J., 2004, MNRAS, 353, 624 
Diemand J., Moore B., Stadel J., 2004, MNRAS, 352, 535 
Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004b, 

MNRAS, 355, 819 
Gao L., De Lucia G., White S. D. M., Jenkins A., 2004a, MNRAS, 

352, LI 

Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel 

J., 1998, MNRAS, 300, 146 
Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel 

J., 2000, ApJ, 544, 616 

© 2005 RAS, MNRAS 000, 000-000 



Gncdin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, 

ApJ, 616, 16 
Governato, F., et al. 2004, ApJ, 607, 688 
Haardt F., Madau P., 1996, ApJ, 461, 20 
Hayashi E., et al., 2004, MNRAS, 355, 794 

Jackson, N., Nair, S., Browne, I. W. A., Wilkinson, P. N., Muxlow, 
T. W. B., de Bruyn, A. G., Koopmans, L., Bremer, M., et al., 
1998, MNRAS, 296, 483 

Katz, N. 1992, ApJ, 391, 502 

Katz N., White S. D. M., 1993, ApJ, 412, 455 

Keeton, C. R., Gaudi, B. S., & Petters, A. O. 2003, ApJ, 598, 138 

Keeton, C. 2002, preprint, astro-ph/01 12350 

Keeton, C. R. 2003, ApJ, 584, 664 

Klypin A., Kravtsov A. V., Valcnzucla O., Prada F., 1999, ApJ, 
522, 82 

Kochanek, C. S. & Dalai, N. 2004, ApJ, 610, 69 

Koopmans, L. V. E., Biggs, A., Blandford, R. D., Browne, 

I. W. A., Jackson, N. J., Mao, S., Wilkinson, P. N., de Bruyn, 

A. G., et al., 2003, ApJ, 595, 712 
Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, ApJS, 111, 

73 

Kravtsov A. V., Gnedin O. Y., Klypin A. A., 2004, ApJ, 609, 482 
Locb A., Peebles P. J. E., 2003, ApJ, 589, 29 
Maccio A. V., 2005, MNRAS, 361, 1250 

Maccio A. V. & Miranda M., 2005, preprint astro-ph/0509598 
Madau P. & Rees M. J., 2001, ApJ, 551, L27 
Mao S., 1992, ApJ, 389, 63 
Mao S., Schneider P., 1998, MNRAS, 295, 587 
Mao S., Jing Y., Ostriker J. P., Weller J., 2004, ApJ, 604, L5 
Mateo M. L., 1998, ARA&A, 36, 435 
Metcalf, R. B. 2002, ApJ, 580, 696 
Metcalf, R. B. & Madau, P. 2001, ApJ, 563, 9 
Metcalf, R. B. & Zhao, H. 2002, ApJL, 567, L5 
Metcalf, R. B. 2004, preprint astro-ph/0412538 
Moore B., Katz N., Lake G., 1996, ApJ, 457, 455 
Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, ApJ, 
499, L5 

Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel 

J., Tozzi P., 1999, ApJ, 524, L19 
Moore, B. 2001, AIP Conf. Proc. 586: 20th Texas Symposium on 

relativistic astrophysics, 586, 73 
Nagai D., Kravtsov A. V., 2005, ApJ, 618, 557 
Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563 
Patnaik, A. R. & Narasimha, D. 2001, MNRAS, 326, 1403 
Reed D., Governato F., Quinn T., Gardner J., Stadel J., Lake G., 

2005, MNRAS, 380 
Ricotti M., Gncdin N. Y., Shull J. M., 2002, ApJ, 575, 49 
Schneider, P. & Weiss, A. 1992, A&A, 260, 1 

Sluse, D., Surdcj, J., Claeskens, J.-F., Hutsemekers, D., Jean, C, 
Courbin, F., Nakos, T., Billeres, M., et al., 2003, A&A, 406, 
L43 

Somcrvillc R. S., 2002, ApJ, 572, L23 
Stadel J., 2001, PhD thesis, U. Washington 

Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New Astronomy, 
9, 137 

Zakharov, A. F. 1995, A&A, 293, 1 

Zentner A. R., Bullock J. S., 2003, ApJ, 598, 49 

Zentner, A. R., Kravtsov, A. V., Gnedin, O. Y., & Klypin, A. A. 

2005a, ApJ, 629, 219 
Zentner, A. R., Berlind, A. A., Bullock, J. S., Kravtsov, A. V., & 

Wechsler, R. H. 2005b, ApJ, 624, 505 



